/**
 * @file   example1.cpp
 * @author StudyBo <studybo@ubuntu18-04>
 * @date   Wed Feb 17 11:08:44 2021
 * 
 * @brief  实现《有限元理论》课程第一章例1的数值计算。
 * 
 * 
 */

#include "Element.h"
#include <Eigen/Sparse>
#include <fstream>

#define pi 4.0*atan(1.0)
typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::Triplet<double> Tri;
typedef Eigen::VectorXd VectorXd;

// g++ -o main example1.cpp -std=c++11 -I /usr/include/eigen3/ -g

double f(double x, double y)
{
    return sin(x);
}

/** 
 * @brief 计算局部系数矩阵的元素值
 * @param ele 当前单元的指针
 * @param i 第一个自由度的局部编号
 * @param j 第二个自由度的局部编号
 * 
 * @return 返回局部系数矩阵的元素 C_ij
 */
double LocalMat(Element<2>* ele, int i, int j)
{
    double ans = 0.0;
    for (int k = 0; k < ele->n_GaussPnt(); k++)
    {
	double xi = ele->GaussionPoint(k)[0];
	double eta = ele->GaussionPoint(k)[1];
	ans += ele->det_Jacobi(xi, eta) * ele->GaussionWeight(k) * ele->phi(xi, eta, i) * ele->phi(xi, eta, j);
    }
    return ans;
}

int main(int argc, char* argv[])
{
    std::vector<double> Nodes({0, pi/4, pi/3, pi/2, 2*pi/3, pi});
    Interval_1_Element* RefElement = new Interval_1_Element();

    // 使用 Eigen库中的稀疏矩阵来保存局部系数矩阵
    SpMat K(Nodes.size(), Nodes.size());
    VectorXd r = Eigen::MatrixXd::Zero(Nodes.size(), 1);

    // 储存所有自由度
    std::vector<Dofs<2> > AllDofs(Nodes.size());
    for (int i = 0; i < Nodes.size(); i++)
    {
	Dofs<2> tmp({Nodes[i], 0.0});
	tmp.SetGlobalIndex(i);
	AllDofs[i] = tmp;
    }
    int n_element = Nodes.size()-1;
    // 局部单元的自由度个数
    int n_Dofs = RefElement->n_Dofs();
    std::vector<Tri> TriList(4*n_element);
    std::vector<Tri>::iterator it = TriList.begin();
    // 按照单元遍历
    for (int k = 0; k < n_element; k++)
    {
	std::vector<Dofs<2> > tmpDofsList(2);
	tmpDofsList[0] = AllDofs[k];
	tmpDofsList[1] = AllDofs[k+1];
	RefElement->SetDofsList(tmpDofsList);
	for (int i = 1; i <= n_Dofs; i++)
	    for (int j = 1; j <= n_Dofs; j++)
	    {
		*it = Tri(RefElement->NdIdx(i), RefElement->NdIdx(j), LocalMat(RefElement, i, j));
		it++;
	    }

	// 计算局部右端项
	for (int i = 1; i <= n_Dofs; i++)
	{
	    double a = 0.0;
	    for (int j = 0; j < RefElement->n_GaussPnt(); j++)
	    {
		double xi = RefElement->GaussionPoint(j)[0];
		double eta = RefElement->GaussionPoint(j)[1];
		a += RefElement->det_Jacobi(xi, eta) * RefElement->GaussionWeight(j) * RefElement->phi(xi, eta, i) * f(RefElement->Global_x(xi, eta), RefElement->Global_y(xi, eta));
	    }
	    r[RefElement->NdIdx(i)] += a;
	}
    }
    // 系数矩阵生成
    K.setFromTriplets(TriList.begin(), TriList.end());
    K.makeCompressed();

    // 使用 Eigen库的共轭梯度求解器求解方程组
    Eigen::ConjugateGradient<Eigen::SparseMatrix<double> > solver;
    Eigen::VectorXd solution;
    solver.setTolerance(1e-12);
    solver.compute(K);
    solution = solver.solve(r);

    // 输出.m文件，可在matlab中画图演示
    std::ofstream os;
    os.open("example1.m");
    os << "x = 0 : pi/25 : pi;\n";
    os << "y = sin(x);\n";
    os << "plot(x, y, 'r--');\n";
    os << "hold on;\n";
    os << "x1 = [";
    for (int i = 0; i < Nodes.size(); i++)
	os << Nodes[i] << ";";
    os << "];\n";
    os << "y1 = [" << solution << "];\n";
    os << "plot(x1, y1, 'b-o', 'MarkerFaceColor','b');\n";
    os.close();
    
}
